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Abstract 

Recent observations by the H.E.S.S. collaboration of the Galactic Centre region have revealed what appears to be 
£3 7-ray emission from the decay of pions produced by interactions of recently accelerated cosmic rays with local molecular 
hydrogen clouds. Synthesizing a 3-D hydrogen cloud map from the available data and assuming a diffusion coefficient 
^1 t)f the form k(E) = kq(E / 'Eq) 5 , we performed Monte Carlo simulations of cosmic ray diffusion for various propagation 
O times and values of and S. By fitting the model 7-ray spectra to the observed one we were able to infer the value 
^ of the diffusion coefficient in that environment (k, = 3.0 ± ^.2kp(?Myr~ x for E = 10 12 5 eV and for total propagation 
M time 10 4 yr) as well as the source spectrum (2.1 < 7 < 2.3). Also, we found that proton losses can be substantial, which 
{S| justifies our approach to the problem. 



1. Introduction 



The High Energy Stereoscopic System (H.E.S.S.) 7-ray 
telescope began operations in 2004 and among its first tar- 
52 gets was the Galactic Centre region [ij. Images obtained 
1 1 from the H.E.S.S. observations had an angular resolution 
0.1°, giving us a map of 7-ray emission of unprecedented 
r ~ H detail. A few point sources were first detected, which 



were compatible with the positions of known objects, like 



o 

> 



^ the black hole Sagittarius A*, the supernova remnant Sgr 
A East and the supernova remnant /pulsar wind nebula 
GO. 9+0.1. By subtracting those known sources, the H.E.S.S. 
^ j collaboration was able to produce a mapping of fainter 
t— H 7-ray emission, stretching over an area of approximately 
OO 300pc by lOOpc around the Galactic Centre [1]. This emis- 
sion seems to follow the contours of molecular gas density, 
as measured by its carbon monosulfide distribution [161 ]. 
up to a distance of approximately 1.3° in galactic longi- 
tude from the Galactic Centre. Beyond that distance the 
7-ray emission diminishes to background levels. 

This apparent emission correlation points towards pos- 
sible models, for which two theories were initially proposed 
[l[. The first postulated that a population of electron 
accelerators produces the observed emission via inverse 
Compton scattering. The objects that would make up such 
a population, such as pulsar wind nebulae, would thrive 
in regions of high-density molecular gas; however the large 
number of sources needed to reproduce the observed emis- 
sion renders this possibility rather unlikely [l|. 

The other theory claims that the 7-rays are produced 
by neutral pion decay, resulting from the interaction of lo- 
cally accelerated cosmic rays with the ambient molecular 
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gas. The lower energy threshold of the H.E.S.S. survey was 
380 GeV, so cosmic ray protons of higher energies would 
be needed to produce the observed 7-rays. Furthermore, 
these cosmic rays would have to have been accelerated 
near the Galactic Centre at some point in the past, yet 
not diffused significantly beyond 1.3° from it. Assuming 
the validity of this theory, one may use this data to in- 
fer the diffusion properties of cosmic rays in the Galactic 
Centre region for various possible sources of cosmic rays. 
Considering as a source of the cosmic rays either the su- 
pernova remnant Sgr A East, with an estimated age of 10 
kyr (l7j . or the black hole Sgr A* with a more remote age 
of activity, Aharonian et al. [1] inferred that the diffu- 
sion coefficient should be no more than Z.hkp(?Myr~ x in 
that area. Busching et al. [3] went one step further with 
an analytical reproduction of the H.E.S.S. observations by 
calculating the emission for different diffusion coefficients. 
Using relativistic protons of mean energy ~ 3TeV to rep- 
resent the cosmic rays and a small number of Gaussian 
functions to represent the molecular clouds they arrived 
at a diffusion coefficient of n = \.Zkp(?Myr~ x . Busching 
and de Jager [4] subsequently expanded their results for 
different ages and on-times for the source. More recently, 
Wommer at al. [18[ employed a different approach, in 
which the motion of individual test particles is computed 
by solving the Lorentz force equation for short time peri- 
ods, and the resulting distributions provide the coefficients 
for the diffusion equation. From their initial conditions for 
the turbulent magnetic field, it is derived that no single 
source can account for the observed emission, but that a 
more continuous source of cosmic ray protons produces a 
better correlation. 

In this paper we present a series of time-dependent sim- 
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ulations of proton propagation in the Galactic Centre in a 
detailed 3D distribution of molecular hydrogen concentra- 
tions for different diffusion parameters. By utilizing a wide 
array of data on molecular clouds in the Galactic Centre 
from various observations we construct a fairly accurate 
density grid of the area of diffusion, as shown in §2.1. We 
can then inject cosmic rays, assumed to be protons, from 
an origin of our choice, which propagate according to the 
diffusion model and its free parameters described in §2.2. 
In §2.3 we deal with the interactions of these protons with 
hydrogen molecules and with the production of photons 
from the resulting neutral pion decays. In §2.4 we briefly 
discuss the likely sources of cosmic ray origin, and the 
implications on their age. In §2.5 we describe the more 
technical aspects of the simulation program, such as the 
initial choice of free parameters and the methods used to 
raise the efficiency of the numerical code. Finally, in §3 we 
present the results of the simulations and in §4 we sum- 
marize and discuss the effectiveness of our approach and 
the likely significance of our results. The present paper 
expands upon the initial results of Dimitrakoudis et al. 

!• 

2. Simulations 

2.1. Synthesizing a hydrogen cloud map 

The area of the Galactic Centre in which we simulated 
the diffusion of cosmic rays is rich in H2 gas, contained in 
a complex setup of high-density clouds, ridges and streams 
that comprise about 10% of our galaxys interstellar molec- 
ular gas, i.e. about 2 to 5 xlO 7 solar masses [16j. Due to 
the high densities involved, tracer molecules were used to 
determine the mass of each gas cloud. To create a realistic 
3-D map of that environment we first obtained the data 
for the 159 distinct molecular cloud clumps, as compiled 
by Miyazaki & Tsuboi [12j], i.e. we assigned to each one 
a galactic longitude, latitude, radius and density. Their 
radial distances are unknown, so we assumed a random 
function to simulate them, within reasonable limits. To 
these data we added the locations and densities of the 
radio-sources Sgr A Q, Sgr Bl, Sgr B2 and Sgr C @, 
which are rich in atomic hydrogen. We then broke down 
the entire CS map of the Galactic Centre by Tsuboi et al. 
[T^ blocks of uniform density, which we treated as 
clouds of equal radius and of random radial distances. Fi- 
nally we added the larger CO clouds by Oka et al. jl3[ in 
order to expand the spatial distribution of our map. Every 
time we added a new set of clouds we checked the possi- 
bility that clouds from previous sets were sharing their 
projected locations with the new ones, and we subtracted 
their masses accordingly. The result was 584 clouds of 
hydrogen gas with a high uncertainty as to their radial 
distances. These clouds were then turned into a 3-D grid 
of 120 x 60 x 60 boxes of uniform density, which form the 
volume of our diffusion model, 5.4 • 10 7 pc 3 . The projection 
of that grid constitutes an area far greater than the extent 



of the observed -ray emission by [1], comprising a total 
mass of 1.2 • 1O 8 M , so particle interactions with clouds 
well outside the observed area of emission are accounted 
for. 

2.2. Diffusion model 

We used the diffusive model of CR propagation, as- 
suming a diffusion coefficient of the form: 

k = k (E/E ) 5 

where E is the energy of the CR protons, Eq = lGeV, 
while fto and S are free parameters, whose values we will try 
to infer through our simulations. The index S is a measure 
of the turbulence of the magnetic field and is assumed to 
be 0.3 < S < 0.6 [15[, while ^0 is a constant whose value 
we will attempt to estimate once 5 has been obtained. The 
diffusion coefficient in our simulations is represented by the 
mean free path i, which is given by the usual expression 

i = 3k/c 

where c is the velocity of light. 

Thus, the test protons move in straight lines of length 
equal to I, after which their directions change randomly. 
At the end of each such free walk, the box number of the 
hydrogen density grid is calculated and a check is made 
for collisions with hydrogen protons. If such a collision 
occurs, the diffusion continues with a lower proton energy, 
as shown in the next section. 

2.3. Production of 7 -rays 

We have assumed that all 7-rays forming the observed 
emission are produced from neutral pion decay. Moreover, 
we have assumed that pions are produced in proton-proton 
collisions through two main channels 

a) p + p — > p + p + 7T + + 7T~ + 7T° 



b) p + p- 



p- 



In case (a) the initial proton loses part of its energy and 
a multiplicity of pions is created that take equal amount 
of the energy lost from the proton. The positive and neg- 
ative pions break down into muons, and ultimately into 
electrons, positrons and neutrinos, while the neutral pions 
decay into 7- rays. Assuming an inelasticity k pp — 0.45 
111 ]. 15% of the initial proton energy will go to the pro- 
duced 7-rays, while the original proton will continue its 
diffusion with 55% of its initial energy. This approach, 
while simplified, gives results which are in good agreement 
with the more detailed spectra produced by Biisching et 
al. 0. 

Channel (b) produces no photons and thus does not 
contribute directly to our observed emission. However, the 
initial proton is turned into a neutron, having lost some of 
its energy in the collision, and will thus continue its propa- 
gation in a straight line, unaffected by the magnetic fields 
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that regulated its trajectory as a proton. It will revert, 
though, back into a proton after half life r = 886.7±0.8«sec 
[T^ . and then it will continue its diffusion as before. Ac- 
counting for time dilation, the distances traveled by neu- 
trons are always much smaller than the mean free path 
for each proton energy (0.028pc for E = 10 12 5 eV; 0.92pc 
for E = 10 14 eV; the radii of molecular gas clouds range 
between 2pc and 40pc), so we can safely assume that this 
case will not affect the diffusion process. 

Using the above we calculate the optical depth for each 
random walk. This is calculated as the product of the in- 
teraction cross section with the mean free path £ and the 
density n of hydrogen molecules. That density is retrieved 
from the grid box where each step ends, a fairly effec- 
tive approximation as most random walks are contained 
within single grid boxes. We have assumed that the cross 
section is a pp = 4 • 10 -26 cra 2 [2], while an increase by a 
factor of 1.3 is sufficient to account for the known chemical 
composition of the Interstellar Medium [10]. The reaction 
probability would then be P = 1 — e _r , where r is the 
optical depth. In practice though, the optical depth is 
a very good approximation of P for low probability val- 
ues (P < 0.15). Thus, to speed up our simulations, we 
can simply define the reaction probability as the optical 
depth, having first checked that its value is always low 
enough to warrant the approximation. At every step, that 
probability is checked against a random number. If the 
random number is smaller than the reaction probability, 
then there is a collision and 7-rays are produced, and the 
proton continues its propagation with diminished energy. 
In the case where its final energy is so low that we are 
no longer interested in its resulting photons, the proton is 
removed from the simulation. 

Apart from proton-proton collisions, CR protons could 
lose energy due to ionization losses, coulomb collisions, 
Compton scattering, synchrotron emission, Bethe-Heitler 
pair production and photo-pion production. Wommer et 
al. have demonstrated that energy loss rates due to the 
four latter processes are insignificant compared to energy 
loss rates from proton-proton collisions [l8[. As for ion- 
ization and Coulomb losses, Mannheim and Schlickeiser 
fl5 | clearly show that their effect is negligible compared 
to that of proton-proton collisions in the local interstellar 
medium. In the higher density environment of the galactic 
centre, and especially within molecular hydrogen clouds, 
their comparative effect would be minimal. 

2.4- Possible sources 

Assuming a single source for the origin of the diffuse 
cosmic rays in the Galactic Centre, then the most likely 
candidates would be either the supernova remnant Sgr A 
East or the black hole Sgr A*. The former has an esti- 
mated age of 10 4 7/r [l7| while the latter could have had a 
burst of activity even further in the past. We have con- 
ducted simulations for both scenarios, selecting an age of 
10 4 yr for the former case and an age of 10 5 yr for the lat- 
ter. In both cases, the galactic coordinates for the source 



are set to I = 0°, b = 0°, which correspond well to the 
location of Sgr A* and to the approximate centre of the 
extended source Sgr A East. 

2.5. Simulation parameters 

Except for the time available for diffusion (which corre- 
sponds to the elapsed time since activity at the source), all 
other parameters are the same for both potential sources. 
We have assumed that the source is active for 100 years 
and that it produces a constant flux of cosmic rays during 
that period. During this time cosmic rays are assumed to 
constantly escape from the source and to start diffusing. 
The proton injection spectrum is divided into six logarith- 
mic bins, ranging from E = 10 12 5 eV to E = 10 15 eV. 
Lower energies would produce 7-rays below the H.E.S.S. 
threshold, while higher energies would have a negligible 
impact on the resulting emission, due to their low numer- 
ical density and their large mean free paths. Each bin has 
the same number of test particles, which is fixed for each 
run, depending on the interaction probabilities for a given 
diffusion coefficient. Those numbers are then weighted af- 
ter each simulation according to a power law distribution 
that fits best the power law of the observed 7-rays from 
H.E.S.S., see [1] (r = 2.29 ± 0.07 s£a£ ± 0.20 sys ). Since the 
weighting procedure is applied to the photon spectrum at 
the end of the simulation, the resulting initial proton spec- 
trum is generated in deference to its modification during 
the diffusion process. 

To improve the efficiency of the simulations, we treat 
the process of continuous production of cosmic rays in the 
following way. Only one burst of cosmic rays of all ener- 
gies is created in the simulation, starting at the beginning 
of the sources activity. The photons produced at the end 
of the simulation are, naturally, the ones observed from 
protons that have propagated for the duration of our sim- 
ulation. In addition to them we also take into account 
all the photons produced in a time period before the end 
of the simulation that is equal to the production time at 
the source (10 2 yr). These correspond to the photons that 
would have been produced at the end of our simulation by 
newer populations of protons produced within that time 
range at their source. The actual production time is of lit- 
tle importance, as long as it much shorter than the prop- 
agation time (as is also seen in [4]). Had we considered all 
the protons to have been produced instantaneously, the 
results would have been the same, provided the number 
of protons was increased to provide the same statistical 
robustness. On the other end, a production time as large 
as 10 3 yr wouldnt significantly alter our results. 

Other parameters besides the two ages are the index 8 
and the proportionality factor ^0 of the diffusion expres- 
sion. For the index S we have chosen the values 0.3, 0.4, 
0.5 and 0.6, while k,o takes on a series of test values by 
increments of log (0.1) around the roughly expected value 
required for a diffusion coefficient that would allow a mean 
propagation of < r > 2 = n- At, in accordance with each to- 
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Figure 1: Reduced x 2 values for different values of «o and S for total 
propagation time 10 4 yr. 

tal time and the mean propagation distance inferred from 
the observations. 

Once each simulation is completed and the results nor- 
malised, the resulting 7-rays are compared against the re- 
sults from H.E.S.S. using the reduced x 2 criterion. 

3. Results 

The resulting reduced x 2 values for cosmic rays origi- 
nating from the supernova remnant Sgr A East, assumed 
to be produced 10 A yr ago, are shown in Fig.l for vari- 
ous values of the diffusion coefficient. The minimum value 
of x 2 calculated is 1.8, and corresponds to S = 0.3 and 
^0 = 0.25/cpc 2 Myr _1 . However we can find minima which 
are less than 2 for all values of 5. If we use those values 
of 5 and ^0 to calculate the diffusion coefficient for pro- 
tons of energy 10 12 5 eV (the lowest energy in our sample 
and also the most important due to their relative abun- 
dance over those at higher energies), we arrive at the re- 
sults illustrated in Fig. 2. We can see that, for each value of 
(5, the diffusion coefficient displays the same minimum at 
k = 3.0 ± Q.2kp(?Myr~ x . This value is higher than that 
calculated by Busching et al. [3[ but close to that sug- 
gested by Aharonian et al. [l[ (less than 3.5fcj9c 2 M?/r _1 ). 

For higher values of k, we also see a disparity between 
the way the reduced y 2 values increase in our results and 
those of Busching et al. The reason for this is that for 
these values the mean free path becomes very large and 
comparable to the dimensions of the area of simulation, 
therefore the simulation becomes inefficient, resulting in 
poor statistics even with increased numbers of test parti- 
cles. 

Fig. 3 shows the reduced x 2 values that correspond to 
a total propagation time of 10 5 yr a possible past activity 
of Sgr A*. If we once again use the values of 5 and ^0 
which correspond to a minimum to calculate the diffusion 
coefficient for protons of energy 10 12 5 eV we arrive at the 



Figure 2: The curves represent the reduced x 2 values for different 
diffusion coefficients, calculated for protons of energy 10 12 - 5 eV, for 
each value of S. We see in all cases a minimum for k « 3kpc 2 Myr -1 . 
The grey line shows the equivalent results from Busching et al. [3] 
for comparison. 

results illustrated in Fig. 4. The resulting best choices for 
the diffusion coefficients appear in Table 1. 
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n[kpc 2 Myr x ] 


0.3 


0.45 


0.4 


0.4 


0.5 


0.56 


0.6 


0.32 



Table 1: Diffusion coefficients for different values of 5 for total prop- 
agation time 10 5 yr. 

This variation for different values of 5 is more pro- 
nounced than in the case of propagation for 10 4 ?/r, but if 
we were to derive a mean value of n like before, we would 
find k = 0.43 ± 0. 05 kp(?Myr~ x . The higher propagation 
time requires the cosmic rays to diffuse more slowly than 
in the case of I0 4 yrs, which is why the resulting diffusion 
coefficient is considerably smaller. 

Fig. 5 and 6 plot the proton spectral indices 7, as these 
were inferred from the simulations, for the two total prop- 
agation times versus ^0 for different values of S. 

In Fig. 5 we can see a fluctuation of 7 values, most 
likely due to irregularities in the interactions of the differ- 
ent energy populations with the complex setup of hydro- 
gen clouds in the small propagation time. A progression to 
lower 7 values as we raise the value of ^0 is, however, evi- 
dent. The values of 7 for the minimum reduced x 2 values 
are in the range 2.1 < 7 < 2.3. For the larger propagation 
time, we can see in Fig. 6 a similar progression to lower 7 
values as a function of k,o. In the range of our best re- 
sults, the proton indices are in the range 2.0 < 7 < 2.1, 
so the proton spectrum is clearly slightly steeper than the 
resulting photon spectrum. These results are in general 
agreement with leaky box model predictions Q. 

Furthermore, it is possible to make an estimate of the 
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Figure 4: Plot of the reduced % 2 as a function of the diffusion co- 
efficient k for different values of S. The total propagation time was 
assumed to be 10 5 yr. Higher values of S tend to favor lower values 
of K. 
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Figure 6: Plot of spectral index 7 inferred from simulations as a 
function of the diffusion normalization for different values of 5. 
The total propagation time was taken to be 10 5 yr. 
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total energy of the protons accelerated at their source, by 
comparing the number of observed photons to the num- 
ber of diffusing protons. Doing so for the lowest observ- 
able proton energy and extrapolating for the energy range 
10 9 - 10 15 eV yields a total energy of (8 ± 1) x 10 49 er#, for 
our different time and optimal diffusion parameters. That 
represents approximately 10% of the energy output of a 
typical supernova explosion. 

The Fermi Gamma-ray Space Telescope should be able 
to observe that area of space in a lower energy range (from 
20 MeV up to 300 GeV). We have repeated our simulations 
for the optimal diffusion coefficients we have found, for a 
total propagation time 10 4 yr, with the proton injection 
spectrum extended to 10 9 eV. In these cases, the emis- 
sion is dominated by the lower energy protons, and thus 
extends from about —0.3° to 0.2° in galactic longitude. 




K[kpc7rvtyr] 



With losses Without losses 



4. Summary /Discussion 

We have presented the results of a series of time-dependent 
simulations of the diffusion of cosmic rays in the Galactic 
Centre region. In the first scenario it was assumed that 
a burst of cosmic rays occurred 10 4 yrs ago, while in the 
second 10 5 yrs ago. Likely candidates could have been the 
SNR Sgr A East or the black hole candidate Sgr A*. Their 
7-ray emission was compared with observations from the 
H.E.S.S. collaboration in order to determine the diffu- 
sion coefficient in that region. For that purpose, a detailed 
3-D map of hydrogen concentrations was synthesized from 
various observations [jj [TH, 0, 3, ■> and two scenarios 
for the origin of the cosmic rays were taken into account. 
For the SNR Sgr A East, 10 4 yrs is the estimated upper 
limit on its age, so the resulting diffusion coefficient should 
be taken as a lower limit. There is much more uncertainty 
concerning the possible age of activity of Sgr A*, but it is 
clear that should that have been in the order of 10 5 yrs ago 
or earlier, the random component of the interstellar mag- 
netic field would have to be considerably more pronounced 
than in the first case. 

The need for such elaborate simulations as were de- 
scribed in the preceding chapters arises from the consid- 
erable losses of protons during their propagation, losses 
which are inextricably connected to the local densities of 
hydrogen gas. In Fig. 7 one can compare the resulting re- 
duced x 2 values from a series of test runs where proton 
energy losses are not taken into account, as opposed to 
the same values from our actual simulations. An underes- 
timation of the diffusion coefficient is evident in the simu- 
lations without losses. Furthermore, this underestimation 
explains the discrepancy between our results and those of 
Biisching et al. [3], as the diffusion coefficient calculated 
without losses is about two times smaller than when losses 
are taken into account in the propagation. 

The above results were derived under the assumption 
that the acceleration of the cosmic rays responsible for 
the observed 7-rays occurred 10 4 or 10 5 yrs ago at a sin- 
gle source, with no subsequent periods of activity at that 



Figure 7: Reduced x 2 values for different values of ^0 and 5 = 0.6 for 
total propagation time 10 4 yr with and without proton energy losses 
taken into account during simulations. 



source. Recent papers [6] have noted that this may be 
a very simplified approach, as there have been many su- 
pernovae in the Galactic Centre region in the past millen- 
nia. Also, the uncertainty in the radial distances of hy- 
drogen concentrations may have had a significant impact 
on the final results. Finally, these simulations assumed 
that the diffusion coefficient remains constant throughout 
the whole region of propagation, and local orderings of 
magnetic fields (including their correlation with molecular 
density within the gas clouds) were not taken into account. 

Those limitations notwithstanding, the results of these 
simulations are useful in providing an estimate of the diffu- 
sion coefficient in the Galactic Centre, taking the different 
magnetic turbulence theories (that become manifest in the 
different values of S) into account. 
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